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I.  INTRODUCTION 


Crustal  inhomogeneities  cause  scattering  of  seismic  waves.  The  effect  is  particularly 
pronounced  for  regional  seismograms  (distance  range  0  -  1000  km)  in  the  frequency  band 
0.5  -  5  Hz.  Surface  wave  phases  such  as  Rg  and  Lg  are  especially  affected;  such  phases 
are  of  potential  importance  for  monitoring  a  Complete  Test  Ban  Treaty  (Blandford,  1981). 
Two  methods  that  can  be  used  to  examine  this  problem  quantitatively  are  physical 
modeling  and  the  use  of  theoretical  calculations  of  the  elastic  wavefield.  In  this  report  we 
attack  the  problem  using  these  methods. 

The  original  thrust  of  the  investigation  was  physical  modeling  at  ultrasonic 
frequencies.  The  models  were  constructed  out  of  nickel  plates,  aluminum  blocks  and  an 
aluminum  powder-epoxy  mixture;  the  metals  have  elastic  velocities  and  densities  similar  to 
igneous  and  metamorphic  rocks  commonly  found  in  the  crust  and  the  aluminum  powder- 
epoxy  mixture  is  an  analog  of  relatively  soft  sedimentary  rocks  such  as  might  be  found  as 
valley  fill  in  the  Basin  and  Range  province.  These  analogies  are  further  strengthened  by 
the  low  attenuation  of  the  metals,  similar  to  igneous  and  metamorphic  rocks,  and  the  high 
attenuation  of  the  aluminum  powder-epoxy  mixture,  similar  to  soft  sedimentary  rocks.  Due 
to  these  similarities,  to  scale  the  model  correctly  it  is  only  neccessary  tc  require  *hat  *he 
ratio  of  the  wavelength  to  the  size  of  the  inhomogeneities  be  the  same  in  the  model  as  in 
the  earth.  Since  wavelength  is  proportional  to  frequency,  if  a  typical  model  freauency  of  1 
MHz  is  to  be  compared  to  a  frequency  of  1  Hz  on  regional  seismograms  then  millimeters  in 
the  model  will  scale  to  kilometers  in  the  earth.  Using  these  principles  models  of  tens  of 
centimeters  across  with  surface  features  a  few  millimeters  in  height  were  constructed.  In 
this  investigation  we  have  concentrated  on  the  scattering  of  Rayleigh  waves  by  surface 
features  such  as  mesas  and  sedimentary  basins;  this  should  provide  an  analog  to  the 
effects  on  the  fundamental  mode  Rg. 


Initially,  under  a  previous  AFGL  contract  (FI  9628-81 -K-0039),  we  attempted  to  build 
complex  models  that  were  analogs  to  real  topographic  and  geologic  situations.  The 
attempt  was  successful,  but  the  complexity  of  the  resulting  seismograms  made  it  difficult 
to  acquire  a  conceptual  understanding  of  the  scattering  process.  Such  an  understanding 
is  important  to  avoid  the  considerable  labor  of  building  a  new  model  for  every  new  area 
examined  and  also  to  allow  estimation  of  effects  where  some  aspect  of  the  geology  such 
as  the  depth  of  a  sedimentary  basin  is  uncertain.  Two  approaches  to  this  problem  were 
explored  in  this  investigation.  One  is  to  examine  simpler  models  in  which  the  scattering 
process  is  less  complicated  and  hence  more  easily  understood.  The  other  is  to  compare 
the  seismograms  from  the  physical  models  with  theoretical  seismograms,  in  this 
investigation  calculated  by  the  finite  difference  technique.  The  theoretical  calculations 
are  inevitably  for  simpler  cases  than  the  experimental  results,  but  nonetheless  can  give 
great  physical  insight  into  the  interpretation  of  the  physical  ultrasonic  modeling. 


II.  ULTRASONIC  EXPERIMENTS 


Introduction 

In  this  chapter  ultrasonic  experiments  carried  out  under  this  contract  and  relevant  to 
It  will  be  discussed.  All  of  the  experiments  examine  Rayleigh  wave  scattering  from  simple 
shapes.  Simple  shapes  are  used  because  of  our  previous  experience  with  more 
complicated  models  where  difficulty  was  encountered  in  interpreting  the  complicated 
seismograms  that  resulted.  Thus  the  focus  of  our  effort  has  been  on  elucidating  the 
underlying  factors  influencing  Rayleigh  wave  scattering  from  topography  by  examining 
simple  shapes  in  conjunction  with  finite  difference  calculations  to  be  presented  in  a 
subsequent  chapter.  First  a  brief  presentation  of  two-dimensional  scattering  from 
rectangular  mesas  is  given;  this  study  was  carried  out  under  a  previous  contract 
(F 1 9628-8 1  -K-0039)  but  is  summarised  here  because  of  its  relevance  to  the  project, 
especially  the  finite  difference  modeling.  Then  some  three-dimensional  models  of  circular 
mesas  and  basins  are  presented. 

Two- dimensional  Mountains 

A  report  has  already  been  given  on  these  experiments,  so  only  a  brief  description  will 
be  given  (see  Nathman,  1980).  The  medium  consisted  of  a  nickel  plate  with  one  edge  cut 
into  the  shape  of  a  rectangular  mesa.  Magnetostrictive  transducers  (Chamuel,  1977)  were 
used  to  couple  into  this  edge;  the  dominant  mode  of  elastic  motion  observed  was  Rayleigh 
waves  with  a  velocity  of  2.776  km/sec  and  a  peak  wavelength  of  nearly  14  mm. 
Receiving  transducers  were  placed  either  between  the  source  transducer  and  the  mesa  to 
receive  the  direct  and  reflected  Rayleigh  waves  or  on  the  other  side  of  the  mesa  from  the 
source  transducer  to  record  the  transmitted  waves  Some  examples  are  shown  in  Figure  i 
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Fig  1:  Rayleigh  waves  reflected  and  transmitted  through  t  wo-din;ensional  mesas. 
Dimensions  are  in  millimeters,  Dir  and  Hef  refer  to  the  p  use  t.-doro  sti -icing  Die 
mesa  and  after  reflection,  observer!  on  the  Same  side  of  the  nes  i  is  th"'  source. 
Trans  refers  to  trie  transmitted  pulse  ouservea  on  the  far  s:-Je  ot  t;:e  mesa  from  the 


source 


The  reflected  and  transmitted  pulse  differ  from  the  source  (  Direct  )  pulse  in  two 
ways.  They  are  of  lesser  amplitude:  Figure  2  shows  the  magnitude  of  the  reflection  and 
transmission  (displacement)  coefficients  as  a  function  of  frequency.  It  should  be  noted 
that  a  considerable  fraction  of  the  incident  Rayleigh  energy  is  scattered  into  body  waves 
at  the  mesa  and  lost  from  the  Rayleigh  wave  field.  The  other  effect  is  the  presence  of 
reverberation,  most  easily  seen  after  the  main  transmitted  pulse.  This  suggests  that  the 
Rayleigh  energy  can  travel  several  different  paths  when  it  interacts  with  the  mesa;  travel 
time  analysis  indicates  that  the  main  transmitted  pulse  seen  in  Figure  1  is  due  to  energy 
that  tunnelled  through  the  base  of  the  mesa,  while  the  later  reverberations  are  due  in 
part  to  energy  that  followed  the  surface  over  the  mesa.  This  question  will  be  examined 
again  in  the  discussion  of  the  finite  difference  results. 

Three- thmensiojuil  Models 

Three-dimensional  experiments  with  circular  mesas  and  basins  have  been  carried  out 
The  experimental  set  up  for  one  of  the  mesa  experiments  is  shown  in  Figure  3.  The  mesa 
was  built  on  the  surface  of  a  circular  aluminum  block  and  was  made  of  an  epoxy-aluminum 
powder  mixture.  The  piezoelectric  source  transducer  produced  a  Rayleigh  wave  pulse  with 
a  dominant  wavelength  of  1.23  mm  (Bullitt  and  Toksoz,  1985)  and  the  receiving  transducer 
was  placed  on  the  other  side  of  the  mesa  so  as  to  provide  a  scattering  angle  (measured  at 
the  center  of  the  mesa)  between  0°  (directly  opposite  the  source)  and  40°  All  of  the 
mesa  and  basin  models  had  the  same  lateral  dimensions  shown  in  Figure  3  with  different 
height  or  depth  of  the  feature.  The  basins  were  either  left  unfilled  or  were  filled  with  the 
epoxy-aluminum  powder  mixture,  simulating  a  low  velocity  sedimentary  fill. 

Seismograms  observed  for  the  model  of  Figure  3  are  shown  in  Figure  4  for  scattering 
angles  of  0°,  10°,  20°.  30°  and  40°  Again  substantial  reverberation  is  observed.  Tra.ol 
times  are  indicated  for  the  following  phases:  rp,  Rayleigh  wave  from  the  source  converted 
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Fig  2  Displacement  reflection  and  transmission  coefficients  as  a  function  of  frequency 
(MHz)  for  case  1 ,  Fig  1 . 


density  2.70  Vp  6.4 


Vs  3.0 


Vr  2.76 


Fig  3:  Geometry  of  three-dimensional  models.  R  indicates  the  position  of  tne  leceiver 
the  source 
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Fig  4:  Seismograms  for  model  of  Fig  3.  Time  is  in  seconds;  output  of  receiving  transducer  is 
shown.  Scattering  angle  is  indicated  on  the  right.  See  text  for  discussion.  Mesa  is 
about  5A  high. 


to  P  at  the  mesa;  rs,  Rayleigh  wave  converted  to  S;  rsr,  Rayleigh  wave  converted  to  S 
upon  striking  the  mesa,  travelling  through  the  mesa  as  S,  converting  back  to  Rayleigh  on 
leaving  the  mesa  (this  wave  type  was  the  dominant  type  for  the  two-dimensional  mesas);  r, 
Rayleigh  wave  that  travelled  over  the  mesa  from  source  to  receiver;  rr,  Rayleigh  wave  that 
reverberated  once  within  the  mesa.  All  of  the  seismograms  are  shown  to  the  same  scale. 
Note  the  loss  of  amplitude  as  the  scattering  angle  increases  and  the  change  in  relative 
importance  of  different  phases  at  different  angles. 

Figures  5-8  show  similar  experiments  for  two  mesas  of  different  heights  and  two 
basins.  In  these  figures  a  control  trace,  i.e.  the  trace  observed  with  no  mesa  or  basin 
present,  is  shown  for  comparison  and  the  traces  are  amplitude  scaled.  For  the  basins,  the 
filled  and  unfilled  (with  epoxy-aluminum  powder  mixture)  are  shown.  In  all  cases 
reverberation  is  observed;  for  the  basins  the  effect  seems  to  be  stronger  for  the  filled 
case  than  the  unfilled  case. 

Summary 

Rayleigh  wave  propagation  across  mesas  and  basins  has  been  observed  for  two-  and 
three-dimensional  ultrasonic  models.  The  height  or  depth  of  the  features  was  of  the  same 
order  as  the  wavelength.  Two  effects  were  consistently  seen.  Firstly,  the  amplitude  was 
diminished  upon  crossing  the  structure.  For  the  two-dimensional  models  it  was  clear  that 
considerable  energy  was  being  scattered  into  body  waves;  presumably,  the  same  was  true 
for  the  three-dimensional  case,  although  signal  to  noise  problems  precluded  the 
observation  of  the  scattered  body  waves.  The  other  important  effect  was  the  presence  of 
reverberations  due  to  converted  and  reflected  waves  within  the  structure. 


sa  is  %  2k  high.  Control  trace  is  the  seismogram  without  the  mesa.  Numbers  next 
ch  trace  are  reduction  factors;  to  find  the  true  relative  amplitude  multiply  the 
amplitude  by  the  factor. 
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EFFECTS  OF  A  1 . 5«m  BASIN,  BOTH 
WITH  CONTROL  TRACE.  AS  A  FUNCTION  OF  ANGIE 


III.  NUMERICAL  MODELING 


fntroductian 

The  finite  difference  method  was  used  to  duplicate  the  experiments  performed  above. 
Inherent  storage  limitations  imposed  by  our  computer  facilities  prevented  us  from  modeling 
three-dimensional  geometries,  but  there  were  aspects  of  the  three-dimensional  models 
which  we  were  able  to  investigate.  The  finite  difference  theory,  applications  and 
references  comprise  a  complete  sub-field  in  applied  mathematics,  and  therefore  will  not  be 
presented  here.  Instead  a  detailed  outline  of  the  relevant  theory,  stability  and  accuracy  is 
given  in  an  Appendix.  We  will  use  this  section  to  discuss  the  results  obtained  through  the 
use  of  the  finite  difference  technique. 

Variations  in  velocity  and  topography  have  significant  effects  on  the  amount  of 
energy  scattered  by  incident  Rayleigh  waves.  In  general,  larger  scale  variations  in 
topography  or  velocity  produce  larger  coefficients  of  reflection  and  smaller  transmission 
coefficients.  Also,  sharp  discontinuities  are  more  efficient  scatterers  of  energy  than 
rounded  features  (Bullitt  and  Toksdz,  1985).  The  Figures  presented  in  this  chapter  are  of 
two  types.  Figure  10,  for  example,  is  a  compilation  of  synthetic  seismograms,  made  by 
storing  both  the  horizontal  and  vertical  nodal  displacements  at  each  time  step.  Figure  1  1 
shows  snapshots  in  time,  where  the  displacements  in  the  grid  are  stored  at  a  single  time 
step.  In  other  words.  Figure  10  represents  the  time  varying  displacement  history  of  a 
series  of  points  on  the  free  surface  and  Figure  1 1  represents  a  spatial  varying 
displacement  map  of  the  entire  model,  at  a  single  instant  in  time. 

There  are  many  subtleties  involved  in  interpreting  both  types  of  diagrams,  in  Figure 
10,  an  outline  of  the  model  being  investigated  is  presented  at  the  bottom  of  the  page  T he 
small  square  boxes  on  the  surface  of  the  model  can  be  thought  of  as  seismometers.  They 


are  the  points  on  the  model  where  the  displacement  histories  were  stored  to  make  the 
synthetic  seismograms.  The  top  part  of  the  diagram  is  a  senes  of  synthetic  seismograms, 
tracing  both  horizontal  and  vertical  components  of  the  displacement  vector.  The  horizontal 
and  vertical  components  are  plotted  at  the  same  scale;  the  horizontal  displacements  are 
smaller  than  the  corresponding  vertical  displacements.  This  difference  is  due  to  the  source 
function,  which  is  prescribed  as  having  approximately  twice  as  much  amplitude  in  the 
vertical  direction  as  the  horizontal  direction.  The  vertical  axis  represents  time,  time  is  zero 
at  the  bottom  and  increases  upwards.  The  horizontal  axis  represents  distance  in  units  of 
wavelength.  Note  the  arrival  time  of  the  incident  Rayleigh  wave,  (R),  migrates  upward  as 
the  Rayleigh  wave  moves  away  from  its  original  location.  Then  after  hitting  the  scatterer, 
it  is  partially  reflected,  (R1).  partially  transmitted,  (R),  and  partially  scattered  to  an  S  wave, 
(s).  The  locus  of  points  connecting  the  first  arrivals  of  a  wave  will  always  be  steeper  for 
slower  moving  waves.  Since  in  a  constant  velocity  medium,  the  incident,  reflected  and 
transmitted  phases  of  a  wave  travel  at  the  same  speed,  the  slopes  of  these  three  arrival 
time  curves  should  all  be  consistent.  This  is  a  useful  criteria  in  identifying  wave  types,  and 
scattering  relationships. 

In  Figure  1 1  the  displacement  fields  are  captured  at  successive  time  steps.  The 
displacement  at  each  grid  point  is  indicated  by  a  line  whose  length  is  proportional  to  the 
displacement  and  whose  orientation  is  in  the  direction  of  the  displacement  vector.  The 
waves  are  marked  using  the  same  convention  described  above.  Note  for  shear  waves  (s), 
where  the  particle  displacement  is  perpendicular  to  the  direction  of  propagation,  the 
horizontal  displacement  amplitudes  are  highest  in  the  part  of  the  wavefront  which  is 
traveling  vertically.  For  P  waves  (p)  the  situation  is  the  opposite,  high  horizontal 
displacement  amplitudes  occur  in  the  part  of  the  wavefront  which  is  propagating 
horizontally.  Rayleigh  waves  (R)  travel  more  slowly  than  P  and  S  waves  and  typically  have 
higher  amplitudes  which  decrease  exponentially  with  depth.  Thus  they  are  most  readily 
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identified  at  or  near  the  surface. 

Step  Discontinuities 

As  mentioned  in  the  previous  chapter,  the  orientation  of  our  investigation  was  to  study 
simple  shapes,  then  after  identifying  the  important  relationships  the  model  was  made  more 
complex.  In  this  vein,  the  step  discontinuity  was  studied  first.  The  coefficients  of 
reflection  and  transmission  were  found  to  depend  on  both  the  ratio  of  the  step  height  to 
the  wavelength  and  the  step  direction.  Six  cases  were  studied,  three  having  a  Rayleigh 
wave  incident  on  a  downstep  and  three  involving  a  Rayleigh  wave  traveling  up  a  step. 
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Downstep: 
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Step  Height  =  (1/2)A: 

Figure  9 
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Step  Height  =  1 X: 

Figures  1 0,  1 1 

■f 

Step  Height  =  (3/2)A: 

Figure  1  2 
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Upstep: 
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Step  Height  =  (1/2)A: 

Figure  1 3 
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Step  Height  -  IX: 

Figures  1  4.  1 5 

Step  Height  =  (3/2)A: 

Figure  16 

In  the  downstep  cases  the  incident  Rayleigh  wave  travels  from  left  to  right.  It  strikes  the 
upper  comer  and  is  partially  transmitted  (R)  partially  reflected  ( R’)  and  partially  scattered 
to  an  S  wave  (s).  The  comer  point  acts  as  a  point  scatterer,  that  is  the  S  wave  moves 
circularly  away  from  the  comer  as  if  the  comer  was  a  source.  The  transmitted  portion  of 
the  Rayleigh  wave  continues  down  the  step  and  scatters  in  the  same  manner  at  the  lower 
comer.  The  resulting  radiation  pattern  of  the  S  waves  is  a  superposition  of  the  waves  from 
the  two  corner  points,  causing  a  complex  interference  pattern.  There  are  lobes  in  the 
wavefront  where  energy  is  concentrated  and  other  portions  of  the  wave  front  where  the 
energy  is  greatly  diminished.  The  location  of  these  lobes  is  a  function  of  step  height: 
Figure  1  1  shows  the  case  of  a  one  wavelength  high  step,  illustrating  the  points  made 


Fig  9:  Synthetic  surface  seismograms  for  a  Rayleigh  Wave  incident  on  a  1/2\ 
homogeneous  downstep  geometry.  Incident  and  transmitted  Rayleigh  wave,  (R). 
reflected  Rayleigh  wave,  (R1),  and  Rayleigh  converted  to  shear  wave,  (s)  are  all 
present. 
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Fig  14;  Synthetic  surface  seismograms  for  a  Rayleigh  Wave  incident  on  a  1  \  homogeneous 
ups t ep  geometry. 


above.  The  amount  of  energy  transmitted  and  reflected  by  the  step  is  directly  related  to 
the  step  height.  As  step  size  increases  less  Rayleigh  wave  energy  is  transmitted,  more  is 
reflected  and  more  is  scattered  into  body  (S)  waves.  These  relations  may  be  seen  by 
comparing  Figures  9,  10  and  12. 

Figures  13-16  show  Rayleigh  waves  climbing  steps  of  various  heights.  Here,  the 
incident  wave  is  traveling  from  right  to  left.  When  the  surface  wave  strikes  the  lower 
edge  of  the  step,  the  nature  of  the  scattering  is  different  from  the  case  of  the  downstep. 
A  larger  proportion  of  the  Rayleigh  wave  is  converted  to  an  S  wave  at  the  first  corner.  This 
S  wave  travels  straight  ahead  into  the  medium  not  evenly  distributed  in  amplitude  about 
the  comer,  but  with  the  majority  of  the  energy  within  r  30°  of  the  horizontal.  The 
transmitted  portion  of  the  Rayleigh  wave  travels  up  the  vertical  face  and  is  again 
scattered  into  an  S  wave  at  the  upper  corner.  Note  in  the  last  frame  of  Figure  1  5  the  S 
wave  (velocity  =  3.2  km/s)  pulls  away  from  the  bottom  of  the  Rayleigh  wave  (velocity  = 
3.0  km/s). 

Comparing  the  two  cases,  the  Rayleigh  wave  reflection  and  transmission  coefficients 
appear  to  be  approximately  equal  for  the  upstep  and  downstep  cases  for  the  same  step 
height.  Therefore,  the  amount  of  energy  scattered  into  the  model  in  the  form  of  body 
waves,  mainly  S,  must  be  equal.  The  radiation  patterns  for  the  two  incidence  directions, 
however,  are  very  different.  Energy  is  Scattered  to  S  waves  by  an  upstep  geometry  in  a 
narrow  band  around  the  horizontal,  while  in  a  downstep  geometry  energy  is  scattered 
approximately  isotropically  from  the  two  corner  points  causing  lobes  of  high  and  low  energy 
around  the  wavefront. 

Infilled  and  hilled  Valleys 

Unfilled  and  filled  valleys  were  investigated  using  both  physical  and  numerical 
techniques.  As  discussed  earlier,  numerical  modeling  was  limited  to  two  dimensions,  while 


physical  modeling  was  extended  to  three  dimensions,  in  comparing  the  results  of  these 
two  modeling  techniques,  it  must  be  kept  in  mind  that  waves  traveling  over  three* 
dimensional  wave  paths  are  not  modeled  by  the  numerical  technique.  However,  many  of  the 
arrivals  observed  in  the  physical  model  are  seen  in  the  numencal  model. 

Valleys  of  four  depths,  (1/2)\,  1  \,  (3/2)X  and  2\,  were  modeled.  An  empty  valley  is 
simply  a  downstep  followed  by  an  upstep.  As  such,  one  would  expect  that  observations 
made  for  simple  steps  would  describe  the  present  problem  and  that  is,  in  fact,  correct.  All 
the  features  described  above  are  evident  in  the  unfilled  valley  geometry,  as  shown  by 
Figures  17  -  21.  The  incident  Rayleigh  wave  strikes  the  downside  of  the  trench,  is 
partially  reflected,  (R1),  partially  transmitted,  (R),  and  partially  scattered  to  an  S  wave  (s) 
The  Rayleigh  wave  continues  this  process  around  all  four  of  the  corners;  each  time  the 
transmitted  Rayleigh  wave  becomes  weaker,  thus  each  scattered  phase  beqdfnes  weaker. 

The  areas  inside  the  black  boxes  in  Figures  22  -  26  were  then  filled  with  low  velocity 
materials  having  S  wave  velocity  0.  -  44  *  0HS,  where  0i<s  is  the  S  wave  velocity  in  the 

half  space.  The  propagation  pattern  of  the  incident  Rayleigh  wave  changes  considerably. 

\ 

The  amplitude  of  the  reflected  Rayleigh  wave  is  much  smaller  if  the  valley  is  filled,  and  the 
conversion  of  surface  wave  energy  to  body  wave  energy  is  more  complex.  There  are  no 
longer  corners  to  act  as  point  scatterers,  and  the  wavefront  entering  the  slow  velocity 
medium  is  severely  distorted.  The  deepest  filled  valley,  Figures  25  and  26,  offers  the 
clearest  view  of  the  scattering  pattern.  There  are  other  important  observations. 

(1)  Note  the  amplitude  of  the  Rayleigh  wave,  (R),  increases  as  the  surface  wave  enters 
the  slower  medium  This  occurs  because  the  wavelength  of  the  surface  wave 
decreases  in  the  slower  medium,  but  the  energy  of  the  wave  stays  approximately 
constant;  this  point  is  made  most  convincingly  in  Figure  26.  This  is  exactly  the  same 
effect  as  the  increase  in  amplitude  of  a  plane  elastic  wave  crossing  an  mterfacp  at 
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17:  Synthetic  surface  seismograms  for  a  Rayleigh  wave  incident  on  a  1  /2K  deep 
unfilled  straight  walled  valley  in  a  homogeneous  medium  Notice  Rayleigh  wave 
reflections  at  the  near,  (R*),  and  far,  (R"),  side  of  the  valley 
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Fig  18:  Synthetic  surface  seismograms  for  a  Rayleigh  wave  incident  on  a  1  X  deep  unfilled 
straight  walled  valley  in  a  homogeneous  medium.  Notice  no  Rayleigh  wave  reflections 
at  the  far  end  of  the  valley  This  is  occurs  because  less  energy  enters  the  basin  due 
to  a  higher  coefficent  of  reflection  at  the  near  side  of  the  valley. 
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Fig  19:  Synthetic  surface  seismograms  for  a  Rayleigh  wave  incident  on  a  3/2\  deep 
unfilled  valley,  in  a  homogeneous  media. 


Fig  21 :  Snapshot  displacements  at  various  times,  for  a  2\  deep  valley,  in  a  homogeneous 
medium. 
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Fig  23:  Same  as  Fig  22,  but  valley  is  1  A  deep. 
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Fig  25.  Same  as  Fig  22,  but  valley  is  2A  deep. 


normal  incidence  from  a  high  impedance  medium  to  a  low  impedance  medium,  well 
known  in  exploration  seismology. 

(2)  Some  of  the  Rayleigh  wave  energy  is  converted  to  S  waves  in  both  the  high  and  low 
velocity  media  [(s)  and  (rs),  Figure  26]. 

(3)  Subsequently,  when  the  body  wave  inside  the  valley  hits  the  lower  boundary,  it  is 
partially  reflected  back  inside  the  valley  (rs'),  and  partially  transmitted  into  the  higher 
velocity  substrate  (rss). 

(4)  From  Figures  22  -  26  it  is  plain  that  a  large  percentage  of  the  energy  in  both  the 
horizontal  and  vertical  components  which  enters  the  filled  valley  is  trapped  and 
reverberates. 

(6)  The  complexity  of  the  wavefield  increases  as  the  valley  becomes  shallower.  It  should 
be  noted  all  the  features  described  above  exist  in  the  shallower  valleys  (Figures  22  - 
24)  but  are  hidden  by  the  multitude  of  different  raypaths. 

Homogeneous  and  fnhomogeneows  Uesas 

Figures  27  -30  show  the  results  for  homogeneous  mesas,  that  is  mesas  made  of  the 
same  material  as  the  substrate.  Similarly  to  the  valley  geometry,  the  mesa  geometry  can 
be  constructed  by  combining  an  upstep  and  a  downstep.  All  of  the  same  scattering 
phenomena  discussed  in  the  step  geometries  are  present  in  the  mesa  geometry.  In  addition 
there  is  a  strong  Rayleigh  wave  produced  when  the  S  wave  scattered  from  the  first  corner, 
which  is  strongly  beamed  in  the  forward  horizontal  direction  due  to  the  upstep,  strikes 
the  opposite  side  of  the  mesa.  This  phase  is  labeled  rsr  in  Figure  29.  The  phase  labelled  R 
travels  around  the  surface  of  the  mesa,  producing  scattered  body  waves  at  each  corner 
The  scattering  pattern  is  more  complex  if  the  mesa  consists  of  low  velocity  material  over  a 
high  velocity  half  space.  Figures  31  -  34,  In  this  case,  the  incident  Rayleigh  wave 
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Fig  27:  Synthetic  surface  seismograms  for  a  1  /2A.  high,  homogeneous,  straight  sided  mesa 
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Fig  32:  Same  as  Fig  3 1 ,  but  mesa  is  1  \  high 
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Fig  34:  Same  as  Fig  31,  but  mesa  is  3/P\  high 


produces  a  larger  reflected  pulse.  In  a  similar  manner  to  the  filled  valley,  the  amplitude  of 
the  transmitted  Rayleigh  wave  increases  as  it  enters  the  low  velocity  material  of  the  mesa. 
There  is  trapping  of  energy  and  reverberation  within  the  mesa,  as  seen  in  Figure  33  This 
occurs  because  most  of  the  energy  which  enters  the  mesa  is  trapped  due  to  the  large 
velocity  discontinuity. 


IV.  COMPARISON  BETWEEN  PHYSICALLY  AND  NUMERICALLY  MODELED  RESULTS 


fntroduction 

The  ultimate  objective  of  this  study  was  to  compare  the  results  obtained  by  the  two 
techniques  presented  in  this  report.  To  achieve  this,  we  plot  a  single  finite  difference 
generated  synthetic  seismogram,  similar  to  those  displayed  in  the  last  chapter,  which  we 
compare  to  the  seismograms  generated  by  the  physical  modeling.  The  two  columns  of 
seismograms  displayed  in  Figure  35  are  seismograms  for  seven  of  the  geometries  we 
investigated.  The  geometries  are,  in  descending  order  from  the  top  of  the  page: 

( 1 )  2A  high  low  velocity  mesa. 

(2)  1  A  high  low  velocity  mesa. 

(3)  The  source  after  traveling  across  a  high  velocity  half  space. 

(4)  Unfilled  1  A  deep  basin. 

(5)  Filled  1  A  deep  basin. 

(6)  Unfilled  2A  deep  basin. 

(7)  Filled  2A  deep  basin. 

In  Figure  35a  on  the  left,  the  traces  show  the  vertical  displacement,  generated  by  the 
finite  difference  method,  at  a  single  seismometer  location  as  a  function  of  time.  Figure  35b 
on  the  right  is  a  compilation  of  the  corresponding  traces  generated  by  physical  modeling. 

In  comparing  the  two  columns  of  traces,  two  traces  generated  by  two  different 
sources  cannot  necessarily  be  directly  compared;  the  seismograms  should  also  be 
compared  with  the  unscattered  source  function  (the  third  trace  in  each  column).  In  both 


columns  of  traces,  the  numbers  which  appear  at  the  beginning  of  some  traces  are 
attenuation  factors.  In  other  words,  a  trace  with  an  attenuation  factor  of  4  has  been 
compressed  by  a  factor  of  four.  Also,  in  comparing  the  two  techniques,  note  that  the 
physical  modeling  was  carried  out  in  three  dimensions,  but  the  numerical  calculations  could 
only  be  done  in  two  dimensions.  The  significance  of  this  is  two-fold:  arrivals  which  come 
from  outside  the  plane  passing  between  the  source  and  receiver  will  not  appear  on  the 
numerically  generated  traces,  and  energy  which,  in  three  dimensions,  leaves  the  plane 
between  the  source  and  receiver  is  not  allowed  to  leave  this  plane  in  the  numerical  models. 

law  Velocity  Mesas 

The  comparison  between  the  two  methods  is  generally  good  for  the  low  velocity  mesa 
geometry.  The  first  four  arrivals  seen  in  Figure  35b  are  present  in  Figure  35a  For  the 
reasons  described  above  the  shapes  of  the  two  seismograms  are  different,  but  the  same 
features  are  present  in  both.  The  early  body  wave  arrivals  are  clearly  seen  on  both 
traces.  Note  that  both  techniques  show  a  broadening  of  the  input  pulse.  This  is  probably 
due  to  high  frequency  components  having  a  lower  transmission  coefficient  for  the  vertical 
step  geometry.  This  trend  is  seen  on  all  the  subsequent  traces.  One  important  difference 
between  the  two  modeling  techniques  is  that  intrinsic  attenuation  is  not  taken  into  account 
in  the  numerical  modeling.  Since  the  low  velocity  composite  we  constructed  is  highly 
absorbent  for  frequencies  in  the  ultrasonic  range,  the  amplitudes  predicted  by  the 
numerical  models  are  consistently  higher  than  those  obtained  from  the  ultrasonic  method. 

Unfilled  and  Filled  VdUeys 

The  lower  four  traces  in  both  columns  were  generated  by  Rayleigh  waves  passing 
through  filled  and  unfilled  valleys.  Both  methods  show  the  delayed  Rayleigh  wave  and  a 
strong  shear  wave  arrival.  It  appears  at  first  sight  that  the  results  from  the  numerical 
technique  are  different  than  those  from  the  physical  modeling,  but  the  differences  are  due 


to  different  plot  scaling,  as  shown  by  careful  measurements.  The  shapes  of  the  resulting 
seismograms  are  different,  but  the  effects  of  the  scatterer  on  the  source  pulse  are 
generally  consistent.  There  is  one  real  difference  between  the  seismograms  generated  by 
a  Rayleigh  wave  which  traveled  over  a  valley  filled  with  low  velocity  sediments.  In  the 
physical  models,  the  peak  amplitude  for  the  filled  valley  is  always  lower  than  the  unfilled 
case.  The  opposite  is  clearly  true  for  the  numerical  models  We  conclude  that  this  is  due 
to  the  lack  of  attenuation  in  the  numerical  modeling  technique. 

Conclusions  and  Predicted  Effects  cm  Regional  Seismograms 

All  the  two-dimensional  scattering  relationships  seen  in  the  physical  modeling  are  also 
present  in  the  numerical  modeling.  Different  models  involving  the  basic  step  geometry  can 
be  investigated  quickly  using  the  finite  difference  technique,  and  this  technique  has  the 
advantage  that  the  displacements,  and  therefore  wavefronts,  can  be  seen  at  all  points  in 
the  media.  The  snapshot  pictures  aid  considerably  in  understanding  the  complex 
scattering  relationships  which  are  present  in  a  simple  two-dimensional  step  geometry. 
Ultrasonic  modeling  has  the  advantage  that  the  geometries  to  be  studied  need  not  be 
restricted  in  size  or  shape  by  available  computer  resources  or  by  limitations  in  the  existing 
finite  difference  theory.  This  technique  also  has  the  advantage  that  attenuation  is 
incorporated  into  the  experiment.  Each  has  important  advantages;  for  this  reason  we  feel 
it  is  important  to  concurrently  pursue  both  types  of  modeling. 

Both  the  ultrasonic  and  the  finite  difference  results  indicate  the  importance  of 
scattering  in  regional  seismic  propagation  and  allow  some  predictions  to  be  marie  about  the 
expected  effects.  For  1  Hz  Rayleigh  waves,  such  as  might  be  found  in  the  important 
regional  phase  Lg,  the  wavelength  is  about  3  km.  Our  results  would  predict  strong 
attenuation  and  reverberation,  leading  to  coda  formation,  from  sedimentary  basins  that 


might  reasonably  be  encountered.  Effects  of  topography  on  such  waves  would  not  be  so 


great  except  in  mountainous  regions,  but  higher  frequency  (shorter  wavelength)  waves 
could  be  strongly  affected  by  normal  topography.  The  studies  especially  emphasise  the 
importance  of  near  surface  low  velocity  materials  on  Rayleigh  wave  scattering  because  of 
the  concentration  of  energy  in  the  low  velocity  material.  The  results  may  be  particularly 
relevant  to  the  Nevada  Test  Site,  which  lies  within  basin  and  range  type  structure.  Finally, 
the  work  reported  here  indicates  why  seismic  station  siting  can  have  such  an  important 
effect  on  the  character  of  regional  events  observed  at  different  stations;  in  the  basin  and 
range,  for  example,  it  would  be  interesting  to  compare  a  station  on  a  mesa  with  one  off  the 


mesa. 


APPENDIX:  NUMERICAL  ANALYSIS 

Finite  Difference:  Theoretical  Development 

Several  numerical  techniques  are  available  for  studying  seismic  wave  propagation. 
The  finite  difference  method  is  one  of  the  most  popular,  as  it  is  relatively  simple  to  program. 
It  is  a  complete  solution  to  the  wave  equation,  so  that  reflected,  refracted,  converted  and 
dispersed  phases  are  included  in  the  solution.  Although  the  finite  difference  method  is 
valid  for  1,2  and  3  dimensions,  we  will  study  two-dimensional  geometries  only,  due  to 
computational  time  and  storage  limitations. 

The  technique  requires  that  displacements  at  all  points  in  the  grid  be  specified  at  two 
successive  time  increments,  separated  by  some  At .  The  displacement  field  is  then 
propagated  through  the  grid  by  successive  discrete  approximations  to  the  wave  equation. 
The  magnitude  of  At  and  the  distance  between  neighboring  grid  points,  Aj  and  A z.  are 
strictly  limited  by  a  stability  equation.  Discussion  of  the  stability  equation  will  be 
postponed  until  a  later  section. 

To  minimize  computational  time  and  storage,  artificial  boundaries  must  be  introduced 
along  the  edges  of  the  grid.  The  correct  specification  of  these  boundaiy  conditions  is 
the  most  difficult  aspect  of  finite  difference  modeling. 

Ihe  himle  Difference  Operator 

A  finite  difference  is  simply  a  numerical  approximation  to  a  continuous  differential 
operator.  We  will  approximate  the  elastic  wave  equation,  but  the  technique  would  be 
similar  for  all  differential  equations  It  should  be  noted  that  the  solution  obtained  by  the 
finite  difference  method  is  not  exact,  there  will  be  rounding  error  due  to  the  discrete 


nature  of  the  computation 
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One  can  use  the  finite  difference  technique  to  obtain  an  approximation  of  the 
derivative  for  any  reasonably  well  defined  function  ?(x).  The  technique  is  to  first  write 
the  function  as  a  Taylor  series.  By  expanding  the  function  in  the  positive  and  negative 
directions  about  some  point  x,  one  can  write; 

?(x  r  Ax)  =  ;r(x)  =  |^-Ax  +  ^-(A x)2r  (Ax)3+0(Ax)4  (1) 

Tnen,  truncating  the  appropriate  series  after  the  second  term  and  rearranging  terms,  the 
first  order  forward  difference  can  be  written  as; 

AxJ-^x}].  (2) 

Similarly,  the  first  order  backward  difference  is  given  by; 

“jJ-  ^  -^-[^(x)-^(x  -Ax)]  (3) 

If  the  expansions  for  the  forward  and  backwards  difference  are  subtracted  before 
truncation,  the  result  is  a  central  difference  formula,  accurate  to  second  order.  By  virtue 
of  sampling  on  both  sides  of  the  point  of  interest,  the  central  difference  formula  is  more 
accurate. 


^(x  +Ax  )-y(x  -Ax )] 


(4) 


The  central  difference  formula  for  the  second  derivative  operator,  d*v°/dz?,  is  also  a 
second  order  approximation.  It  is  most  easily  found  by  first  forward  differencing,  then 
backward  differencing. 


d_ 

dx 


1 £l 

dx 


(5) 
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The  truncated  Taylor  series  can  then  be  written  as; 


92-,-0  ^  1 

"  (Ax) 


-^[^(x  -t-Ax)-2v5(x)+^(x  -Ax)] 


(6) 


It  should  be  noted,  as  is  the  case  with  mixed  continuous  derivatives,  the  order  of 
differencing  has  no  effect  on  the  outcome.  Thus,  equation  (6)  could  have  been  obtained 
by  backward  differencing  followed  by  forward  differencing. 

A  more  theoretical  discussion  of  the  ideas  presented  above  is  discussed  in  Burden  el 
al.  (1981).  There  are  many  papers  on  the  application  of  the  finite  difference  method  to 
wave  propagation  problems;  the  work  of  Alterman  and  Karal  (1968),  Boore  (1972),  and 
Kelly  et  ai.  (1976)  is  often  considered  the  most  important  early  work  in  this  discipline. 


I 


Application  of  the  finite  Difference  Method  to  the  Elastic  Wave  Equation 


If  gravity  forces  are  ignored,  the  elastic  wave  equation  for  an  isotropic, 
inhomogeneous  media  is  given  by 

32,7  f  1  f  1 

iZj  -  V  x(MV  xd)  *  2^(VM  'T)i2+(VmK  if  x  ( V  x  u  )j  (7) 

where  p  is  density,  if  is  the  displacement  vector,  \  and  p  are  Lame's  parameters  and  t  is 
time.  Let  i  and  4  be  the  horizontal  and  vertical  directional  unit  vectors  in  rectangular 
coordinates,  with  4  being  positive  down.  For  the  two-dimensional  case  in  a  perfectly 
elastic  and  isotropic  body,  the  wave  equation  can  be  reduced  to  two  coupled  equations  in 
terms  u  and  v,  the  horizontal  and  vertical  components  of  the  displacement  vector  (Ewing 
et  al..  1957): 
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These  equations  fully  describe  the  motion  of  compressional  (P)  and  vertically  polarized 
(SV)  waves  within  the  media.  There  is  no  need  to  consider  the  horizontally  polarized  shear 
(SH)  waves  in  a  two-dimensional  analysis  because  SH  motion  is  orthogonal  to  both  P  and  SV 
motion.  This  orthogonality  prevents  coupling  of  SH  energy  to  SV  and  P  waves. 


The  Homogeneous  Case 

If  the  physical  properties  of  a  medium  are  invariant,  equation  (8)  can  be  simplified  and 


written  as 


d2! 1 


dt 


=  a2 


d2^  2 

p — tr  -  a 
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d*~U  rj^w 
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d^u 

d^u. 
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d^u 
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dz2 
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(9) 


where  a  and  8  are  the  P  and  S  wave  velocities,  a  =  ~\/  -+—&-  and  B  =  -\ /  Since 

v  P  v  P 


Lame's  parameters  are  not  differentiated  at  each  node,  this  form  of  the  wave  equation  is 
much  faster  to  compute. 


For  finite  difference  calculations,  the  continuous  derivatives  in  equation  (9)  are 
replaced  with  their  discrete  counterparts: 


u(z  ,x,t  +  At )  =2(1  -a2-^2)  u(z  ,x,t )  -  u(z  ,x ,f  -AO  +  At  ^ 

x  [a2 (11(2  ,x  + 1  ,t  )ru (z  ,x  -1  ,t )) 

+  /?2(u(z  +1  ,x,t  )+u(z  -1,x,t )) 

2  2 

+■  — ^  -  [n(z  +  1,x  +  1.t)  -  xl(z  +1  ,x -1  ,t )  - 
ix(z-1,x+1,t)  +  tx(z  -1  ,x  -1  ,t  )]] 


u  (2  ,x,t  +  Af )  =  2  ( 1  -a2-/?2)  u  (z  ,x,t )  -  u,( z ,x,t  -At )  +  At 

x[^2(tL(z,x+1,t)  +  n,(z  ,x  -1  ,t )) 

-I-  o2(u,(z  +1 , x,f )  +  u,(z  -1  ,x,f )) 

2  2 

+•  -  ^  [u(z  +  1  ,x  + 1  ,t )  -  u(  z  + 1  ,x -1  ,t )  - 

u(z  — 1  ,x  + 1  ,t )  +■  11  (z  —  I  ,x -1  ,t )]] 


(10) 


These  equations  constitute  an  explicit  finite  difference  scheme  That  is,  the 
displacements  u  and  u  at  each  point  in  the  gr.d  can  be  calculated  at  time  step  (.*  *■  At) 


V- V 


using  only  the  displacements  at  the  two  previous  time  steps  t  and  (t  -  It). 


The  fnhomogeneaus  Cdse 


For  an  inhomogeneous  medium  Lame's  parameters  and  density  as  well  as  the 
displacements  vary  spatially.  Accordingly,  not  only  displacements  but  also  these  physical 
properties  must  be  differenced  at  ail  points  in  the  grid.  This  adds  significantly  to  the 
computational  burden.  For  equal  sized  grids  CPU  time  requirements  for  the  inhomogeneous 
formulation  were  between  two  and  three  times  greater  than  for  the  homogeneous  case.  To 
reduce  run  time  and  storage  requirements  we  assume  that  density  is  constant  everywhere 
in  the  grid.  By  rearranging  terms,  expressions  involving  \  and  ^  can  be  replaced  by 
expressions  involving  a  and  /?.  Derivatives  involving  these  velocity  variables  are  solved 
using  the  chain  rule.  Thus, 


a  2j  \ 

—  aiz  ,x )  — 
ox  ox 


becomes 


(12) 

dx  dx  Qx  *■ 

Terms  involving  mixed  derivatives  are  solved  in  a  similar  way.  Evaluating  equation  (8)  with 
this  technique,  and  replacing  continuous  derivatives  with  the  corresponding  discrete 
approximations  yields: 


u(z,x,t  +St)  =  2  u(z,z,0  -  u(z,x,t  -St )  4-  M^y. 


x  [  — — ry{(a^(z  ,x  +1)  4-  o3[z,z))(u(z  ,z  4-1 ,0  -  v  (z.z.O)  - 
2Sx^ 

(a^z  ,x  -1 )  +  a^(z  ,x  ))(u(z  ,x,t )  -  u(z  ,x  -1  ,t ))] 

4- — T(a^(z^r+1)  -  2/?^(z  ,x  +1  ))(it,(z  +1  ,x  +1  ,t)  -  xl(z  -1  ,x  +1  ,t  ))- 

4A z  Az 
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4Ax^z 
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4-  — — p{(0^.Z  4-1  ,x)  4-  0^Z  ,z))(u(z  4-1  ,2,0  -  u(z  ,Z,0)  - 
2Sz^ 

(/?3(z-1,z)  4-  p2(z  ,x))(u(z  ,x  ,t )  -  u(z -1  ,z,0)]] 


U,(z  ,2,i  +00  =  2  12  (  2,1,0  -X i,(z, X  A  -AO  4-  St 
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Accuracy.  Boundaries  and  the  Source 


Accuracy 

Inaccuracy  in  finite  difference  calculations  is  primarily  due  to  grid  dispersion  and 
inadequate  stability  at  boundaries.  If  Ax  or  Az  is  too  large,  the  low  frequency  component 
of  compressional  waves  will  travel  faster  across  the  grid  than  the  high  frequency 
components.  This  grid  dispersion  is  purely  a  numerical  artifact  which  can  be  minimized  by 
keeping  grid  spacings  small.  Grid  dispersion  of  shear  waves  is  much  more  complex.  It  has 
been  shown  that  for  some  combinations  of  Poisson's  ratio  and  propagation  direction,  high 
frequencies  will  travel  faster  than  low  frequencies  (Trefethen,  1982;  Bamberger  et  al., 
1980). 

Sato  and  Ishihara  (1983)  studied  the  effects  of  grid  size  on  the  velocity  of  various 
waves  traveling  in  a  finite  difference  grid.  They  found  that  the  high  frequency  components 
in  compressional  waves  travel  more  slowly  than  the  low  frequency  components.  For  shear 
waves,  the  behavior  is  reversed.  In  addition,  they  obsen/ed  that  all  waves  appear  to 
travel  faster  at  angles  near  45°  to  the  grid.  This  is  due  to  the  fact  that  the  waves  are 
sampled  less  frequently  as  the  distance  between  nodes  becomes  greater,  as  happens  at 
45°.  The  amount  of  dispersion  depends  also  on  the  ratio  A J  p;  dispersion  is  minimized  when 
this  ratio  is  unity. 

The  magnitude  of  A t  has  a  significant  effect  on  grid  dispersion  and  the  amount  of  high 
frequency  information  which  is  preserved  This  is  clearly  pointed  out  by  examining  the 
effects  of  discrete  differentiation  on  the  harmonic  function, 

(14) 

Differentiating  twice  yields 
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Using  the  second  order  central  difference  formula  given  in  equation  (6) 


dZAt)°  ^  e  ^  -2e°+exu* 

(T7F 


and  applying  trigonometric  identities  reduces  equation(16)  to 


(16) 


dZAt)° 


?cos:uAt  -2 


(17) 


where  the  apparent  frequency  £>  is  given  by; 


=  2_sin(HMi) 
At  2 


2:jA  t 
2At 


(18) 


This  approximation  holds  when  ^Af  is  small.  It  should  be  noted  that  for  high  frequencies 
the  apparent  frequency  2  becomes  significantly  smaller  than  the  true  value  Thus,  high 
frequency  waves  will  be  dispersed  more  strongly  than  low  frequency  waves  given  the 
same  At-  This  relationship,  shown  in  Figure  A1,  indicates  that  frequencies  greater  than 
1  /  At  will  not  be  adequately  approximated  by  the  finite  difference  technique. 


The  effects  of  the  spatial  grid  size  on  stability  can  be  investigated  by  analyzing  the 
acoustic  wave  equation 


1 

dt'Z  dz^  dy^ 


(19) 


Applying  a  three-dimensional  Fourier  Transform,  we  obtain 


(20) 


Applying  the  second  order  central  difference  approximation  to  equation  ( 1  9), 


a  uAf  .  4  ,  2,  4  a  VV, 

2-  sinn  — — )  -  sinn  — - — )  -  7- — -7  sinn  — n — ). 
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(lif) 

Rearranging  terms  equation  (21 )  can  be  written  as; 


(Ay)2  *  2 


(21) 


axAt.  ?( Af )2  r  2/  \  (Ax)2  a  fcvAy 

smn  — -)  =  v^-—-7  [sinn  — — )+  7— -7  sinn  -*r — )]■ 


(22) 


(Ax)^  2  '  (Ay )' 

Taking  the  maximum  value  of  the  sine  as  1  and  setting  Ax  =Az ,  equation  (22)  can  be 
rewritten  as 


Af  < 


1  Ax 
x/5  v 


(23) 


For  the  elastic  case,  this  relationship  has  to  be  modified  slightly  to  include  the  shear  wave 
velocity. 


Af  < 


min(Ax,Az ) 


(24) 


Vn2*/?2 

Equations  (23)  and  (24)  are  the  stability  conditions  for  two-dimensional  acoustic  and 
elastic  finite  difference  schemes.  In  equation  (24),  a  and  p  are  taken  to  be  the  largest  P 
and  S  wave  velocities  in  the  grid,  Ax  and  Az  are  the  distances  between  grid  points  in  the 
x  and  z  directions  and  Af  is  the  time  step  necessary  to  calculate  the  time  derivative 
accurately.  Since  the  shear  wave  velocity  is  always  assumed  to  be  less  than  the  P  wave 
velocity,  one  can  see  that  the  elastic  form  of  the  stability  equation  is  somewhat  more  lax 
than  the  acoustic  form.  These  stability  equations  have  been  shown  to  be  adequate  for 
most  explicit  finite  difference  approximations  of  both  the  acoustic  and  the  elastic  wave 
equations  (Kelly  et  al..  1  9  76). 

From  the  relationship  between  frequency,  velocity  and  wavelength  the  maximum 


allowable  grid  size  in  a  finite  difference  mesh  is  given  by 


(25) 


Ax 


nf  max 


where  n  is  the  number  of  points  per  wavelength,  /nwis  the  maximum  frequency,  and  i'rnn 
is  the  minimum  velocity  value  in  the  model  (the  slowest  Rayleigh  wave  velocity  occurring  in 
the  model). 


One  would  like  to  maximize  the  number  of  points  per  wavelength  in  order  to  preserve 
as  much  high  frequency  information  as  possible.  Doing  so,  however,  forces  At  to  be  very 
small,  thus  increasing  the  number  of  time  steps  necessary  to  propagate  a  wave  a  given 
distance.  This  in  turn  increases  computational  time  and  expense.  Our  analysis  indicated 
that  20  points  per  wavelength  was  sufficient  to  maintain  adequate  high  frequency 
information  for  our  purposes.  Choosing  20  points  per  wavelength  and  Ax  =  Az  =  At  =  1.0 
required  that  a<  .7,  £<a/  s/3. 

Free  Surface  Boundaries 

Our  investigation  studied  the  effects  of  crustal  inhomogeneities  on  the  passage  of 
surface  waves.  Therefore,  a  free  surface  boundary  at  the  top  of  our  grid  is  needed.  By 
definition,  the  normal  and  tangential  stresses  must  vanish  at  a  free  surface  (Ewing  et  ai, 
1  957).  For  a  horizontal  free  surface,  this  is  represented  by 


(a2_2i9z>au 


4-  'X2^-  =  0 
dx  dz 


(26) 


and 


du 

dz 


(27) 


To  approximate  the  horizontal  free  surface,  an  additional  row  of  grid  points 
(pseudonodes)  is  introduced  above  the  free  surface.  Displacements  at  the  pseudonodes 
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are  found  using  forward  differences,  while  the  displacements  at  the  interior  rows  are 
determined  with  central  differences.  This  will  increase  accuracy  at  the  nodes  within  the 
media.  The  same  general  technique  is  used  for  vertical  free  surfaces,  except  equation 
(26)  is  rotated  by  rr/  2. 

If  the  free  surface  contains  90°  or  270°  comers  the  displacements  at  the  comer 
points  are  found  using  the  central  difference  formulation  and  the  displacements  at  the 
points  adjacent  to  the  comer  points  are  determined  using  a  double  pseudonode.  This 
technique  is  described  in  more  detail  by  Fuyuki  and  Matsumoto  (1980). 

Other  Boundaries 

Time  spent  computing  a  finite  difference  grid  increases  rapidly  as  additional  nodes  are 
appended  to  the  grid.  For  this  reason,  it  is  desirable  to  make  the  grid  as  small  as  possible, 
without  interfering  with  the  topology  under  study.  To  accomplish  this  objective,  boundaries 
are  needed  at  the  edges  of  the  grid  which  appear  transparent  to  outgoing  waves,  and 
therefore  prohibit  inward  reflections. 

There  are  three  common  techniques  which  have  been  invoked  to  achieve  the  desired 
behavior  at  the  boundaries.  The  first  technique  (Lysmer  and  Kuhlemeyer,  1969)  adds  a 
viscous  damping  term  to  the  wave  equation,  to  damp  the  normal  and  shear  components  of 
the  stress  tensor  along  the  boundaries.  This  technique  works  well  for  compressional 
waves  at  all  angles,  but  does  not  completely  attenuate  the  reflected  shear  waves  An 
alternative  method  was  put  forth  by  Smith  (1974),  in  which  the  full  wave  equation  was 
solved  twice,  once  with  Dirichlet  boundary  conditions  and  once  with  Neumann  boundary 
conditions.  The  two  solutions  are  then  summed  to  completely  annihilate  reflections  from  all 
incidence  angles.  The  basis  of  the  technique  is  the  result  that  waves  reflecting  from  zero 
displacement  and  zero  stress  boundaries  are  equal  in  magnitude  but  opposite  in  sign. 
Clearly  this  technique  has  exactly  the  attributes  that  are  desired,  but  the  computational 
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burden  of  solving  the  full  wave  equation  twice  at  each  boundary  is  too  great.  Clayton  and 
Engquist  (1977)  derived  a  formulation  which  does  not  completely  eliminate  reflections  from 
all  angles,  but  takes  only  the  same  number  of  operations  as  a  point  on  the  interior  of  the 
grid.  They  used  the  paraxial  approximation  to  the  full  wave  equation  to  separate  the  in¬ 
going  and  outgoing  components  of  the  wave.  By  doing  so,  the  in-going  component  can  be 
reduced  at  the  boundary.  The  combination  of  efficiency  and  adequate  accuracy  led  us  to 
chose  the  Clayton-Engquist  formulation  for  all  of  our  numerical  experiments. 

The  Source 

Historically,  the  Ricker  wavelet  (Ricker,  1977)  has  been  the  most  common  source 
function  for  numerical  modeling  of  surface  waves  because  the  displacement  field  is 
localized  in  both  the  space  and  wavenumber  domains.  In  the  space  domain  the 
displacement  at  the  free  surface  is  depicted  by  Figure  A2.  In  the  wavenumber  domain  the 
spectrum  for  the  Ricker  wavelet  is  given  by 

S(k)  =  U/fca)3axp[1  -Ob/fc,)2]  (28) 

where  k  is  wavenumber  and  k0  is  the  maximum  amplitude  wavenumber.  Equation  (28)  is 
shown  graphically  in  Figure  A3. 

In  this  study  we  define  the  displacement  at  time  zero  at  the  free  surface  as  the 
Ricker  wavelet  and  calculate  the  displacements  at  all  points  below  the  surface  using  the 
eigenfunction  for  a  Rayleigh  wave  on  a  half  space.  The  source  calculation  was  performed 
in  the  wavenumber  domain  and  transformed  to  the  space  domain  by  Fourier  transformation. 
Thus,  the  initial  displacements  for  all  points  in  the  grid  can  be  calculated  from 

il(z,x,n=  jrJ'E(k)S(k)dk  (29) 

Here  EX k )  is  the  half  space  eigenfunction  for  Rayleigh  wave  displacement 


where 


Elk)  = 
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ik 

-kRc 


kIR, 

ik/ 


exp(  -A:/4  z ) 
exp(  -kRsz) 


exp(i k{x  -  I'.f )) 


Rf=  1 


/  is  the  forward  wave  impedance  and  V  is  the  half  space  Rayleigh  wave  velocity 
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